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EXECUTIVE SUMMARY 


The U.S. Department of Energy (DOE) and the Nevada Division of Environmental 
Protection (NDEP) have reached agreement on a corrective action strategy applicable to address 
the extent and potential impact of radionuclide contamination of groundwater at underground 
nuclear test locations. This strategy is described in detail in the Federal Facility Agreement and 
Consent Order (FFACO, 2000). As part of the corrective action strategy, the nuclear detonations 
that occurred underground were identified as geographically distinct corrective action units 
(CAUs). The strategic objective for each CAU is to estimate over a 1,000-yr time period, with 
uncertainty quantified, the three-dimensional extent of groundwater contamination that would be 
considered unsafe for domestic and municipal use. 

Two types of boundaries (contaminant and compliance) are discussed in the FFACO that 
will map the three-dimensional extent of radionuclide contamination. The contaminant boundary 
will identify the region with 95 percent certainty that contaminants do not exist above a threshold 
value. It will be prepared by the DOE and presented to NDEP. The compliance boundary will be 
produced as a result of negotiation between the DOE and NDEP, and can be coincident with, or 
differ from, the contaminant boundary. 

Two different thresholds are considered for the contaminant boundary. One is based on the 
enforceable National Primary Drinking Water Regulations for radionuclides, which were 
developed as a requirement of the Safe Drinking Water Act. The other is a risk-based threshold 
considering applicable lifetime excess cancer-risk-based criteria 

The contaminant boundary for the Faultless underground nuclear test at the Central Nevada 
Test Area (CNTA) is calculated using a newly developed groundwater flow and radionuclide 
transport model that incorporates aspects of both the original three-dimensional model 
(Pohhnann et al, 1999) and the two-dimensional model developed for the Faultless data decision 
analysis (DDA) (Pohll and Mihevc, 2000). This new model includes the uncertainty in the three- 
dimensional spatial distribution of lithology and hydraulic conductivity from the 1999 model as 
well as the uncertainty in the other flow and transport parameters from the 2000 DDA model. 
Additionally, the new model focuses on a much smaller region than was included in the earlier 
models, that is, the subsurface within the UC-1 land withdrawal area where the 1999 model 
predicted radionuclide transport will occur over the next 1,000 years. 

The purpose of this unclassified document is to present the modifications to the CNTA 
groundwater flow and transport model, to present the methodology used to calculate contaminant 
boundaries, and to present the Safe Drinking Water Act and risk-derived contaminant boundaries 
for the Faultless underground nuclear test CAU. 
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1.0 INTRODUCTION 

The U.S. Department of Energy (DOE) has a responsibility to assess the potential health 
risk to workers and the public from possible exposure to environmental contamination resulting 
from nuclear testing, and the Nevada Division of Environmental Protection (NDEP) is the 
regulatory agency responsible for protecting public health and the environment in Nevada. These 
organizations reached agreement on a corrective action strategy applicable to address the extent 
and potential impact of radionuclide contamination of groundwater at underground nuclear test 
locations. This strategy is described in detail in Revision 1 to Appendix VI, Section 3, of the 
Federal Facility Agreement and Consent Order (FFACO, 2000). This corrective action strategy 
also has undergone an independent technical peer review (ASME/IRS, 2001). 

As part of the corrective action strategy, the nuclear detonations that occurred underground 
were identified as geographically distinct corrective action units (CAUs). The objective for each 
CAU is to estimate over a 1,000-yr time period, with uncertainty quantified, the three- 
dimensional extent of groundwater contamination that would be considered unsafe for domestic 
and municipal use. The quantification of the uncertainty is derived from an uncertainty analysis 
of the flow and transport model. This report identifies various methodologies that could be used 
to calculate this region of contaminated groundwater. 

The contaminant boundaries for the Faultless underground nuclear test at the Central 
Nevada Test Area (CNTA) are calculated using a newly developed groundwater flow and 
radionuclide transport model that incorporates aspects of both the original three-dimensional 
model (Pohhnann et al., 1999) and the two-dimensional model developed for the Faultless data 
decision analysis (DDA) (Pohll and Mihevc, 2000). The original three-dimensional model was 
approved by the State of Nevada conditional on additional uncertainty analysis, which was 
performed in the DDA. To be true to the State's intent for a more rigorous incorporation of 
uncertainty, DOE/Nevada Operations Office decided to enhance the original three-dimensional 
model with the uncertainty features developed in the DDA before proceeding with determining 
contaminant boundaries for Faultless. This new model includes the uncertainty in the three- 
dimensional spatial distribution of lithology and hydraulic conductivity from the 1999 model as 
well as the uncertainty in the other flow and transport parameters from the 2000 DDA model. 
Additionally, the new model focuses on a much smaller region than was included in the earlier 
models; that is, the subsurface within the UC-1 land withdrawal area where the 1999 model 
predicted radionuclide transport will occur over the next 1,000 years. 

Numerical Monte Carlo methods are utilized to incorporate the uncertainty in the stochastic 
groundwater flow and radionuclide transport model. Calibration of the model is addressed 
through application of Bayesian methods to weight each realization based on the degree to which 
the simulated flow field matches the flow field observed at Faultless. 

Two types of boundaries are discussed in the FFACO that will map three-dimensional 
groundwater regions, a contaminant boundary and a compliance boundary. Appendix VI of the 
FFACO (2000) provides the following description of these boundaries: 

“A contaminant boundary is the model-predicted perimeter which defines the extent 
of radionuclide-contaminated groundwater from underground testing above 
background conditions exceeding the Safe Drinking Water Act (SDWA) standards. 

The contaminant boundary will be composed of both a perimeter boundary and a 
lower hydrostratigraphic unit boundary. The computer model predicts the location of 
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this boundary within 1,000 years and must do so at a 95% level of confidence. 

Additional results showing contaminant concentrations and the location of the 
contaminant boundary at selected times will also be presented. These times may 
include the verification period, the end of the five-year proof of concept period, as 
well as other times that are of specific interest. 

From the contaminant boundary predicted by the computer model, a compliance 
boundary will be negotiated between NDEP and DOE. The compliance boundary 
will define the area within which the radiological contaminants above the SDWA 
standards relative to background are to remain. DOE will be responsible for ensuring 
compliance with this boundary. The compliance boundary may or may not coincide 
with the contaminant boundary. If the predicted location of the contaminant 
boundary cannot be accepted as the compliance boundary, an alternate compliance 
boundary will be negotiated by both parties.” 

The water-quality standard used for detennining whether groundwater poses a health risk is 
specified in the FFACO as the National Primary Drinking Water Regulations (NPDWRs) for 
radionuclides, which were developed as a requirement of the Safe Drinking Water Act. Though 
not required by the FFACO, consideration should also be given to the criteria that the U.S. 
Environmental Protection Agency (EPA) considers in establishing regulatory standards for 
drinking water contaminants. The EPA target range for lifetime excess cancer risk is not to 
exceed 10' 4 (1/10,000) and ideally is less than 10" 6 (1/1,000,000) (EPA, 2000a), and 
corresponding maximum contaminant levels (MCLs) for specific radionuclides. 

This report is organized to first present the specifics of the revised groundwater flow model 
and radionuclide transport model. The approach to defining a contaminant boundary is then 
discussed, beginning with the process of determining what constitutes unsafe water, then 
proceeding with the process of incorporating uncertainty in the boundary. Examples are shown 
and discussed. Finally, the contaminant boundary for the Faultless test is presented. 
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2.0 FLOW MODEL 

2.1 Conceptualization of Groundwater Flow 

There have been no additional data collection efforts at Faultless since the completion of 
the 1999 model and the conceptualization of groundwater flow has not changed since that time. 
A brief overview of the conceptual flow model as described by Pohhnann et al. (1999) is 
included here to reacquaint the reader with the principal components. 

The flow model is based on hydrogeologic conditions prior to the Faultless test, under the 
assumption that transport over the long term is controlled by these factors rather than the 
relatively short-tenn effects of the test. Furthermore, flow is considered to be at steady state, 
owing to the large size of the Hot Creek Valley hydrologic system and the absence of excessive 
groundwater withdrawals. Local structural features such as faults are not explicitly included due 
to the lack of information regarding their subsurface locations and hydraulic characteristics. 

The principal hydrogeologic units present in the UC-1 land-withdrawal area are divided 
into three categories, based on their lithologic and hydraulic properties. Category 1 represents 
Quaternary alluvium containing pebble- to boulder-size fragments enclosed in a clay-cemented 
matrix. Category 2 consists of Tertiary tuffaceous sediments, bedded tuffs, and partially welded 
tuffs that exhibit various degrees of alteration (argillation, zeolitization) and thus larger amounts 
of clay matrix materials than the other two categories. As a result, these rocks have the lowest 
hydraulic conductivity of the three categories, as demonstrated by the results of the CNTA 
straddle packer testing (Pohhnann et al., 1999). Tertiary rhyolites and densely welded tuffs 
comprise Category 3. These rocks are highly fractured and where spatially connected are 
considered the primary pathways for groundwater flow and transport at Faultless. 

Hydraulic head measurements in wells in Hot Creek Valley indicate that groundwater flow 
in the alluvium is directed toward the south, generally following the slope of the valley surface. 
Head relationships in the deeper volcanic system, including surrounding valleys, indicate flow 
toward the northeast and east toward regional discharge points further south in Hot Creek Valley 
and east in southern Railroad Valley. Strong vertically downward hydraulic gradients are present 
north of UC-1, where land surface elevations are higher and recharge from precipitation is likely 
to be taking place. To the south, strong vertical hydraulic gradients from the volcanic section 
upward toward the alluvium are present and may be related to the regional discharge. Head 
measurements at HTH-1 prior to Faultless suggest a very low downward-directed ambient 
gradient. 

2.2 Flow Modeling Methods 

2.2.1. Solving the Flow Problem 

The flow modeling closely follows the process established in the 1999 model and though a 
few enhancements were made to the model in terms of approaches, no new data were 
incorporated. As was done for the 1999 model, uncertainty and heterogeneity in hydraulic 
conductivity ( K) are incorporated through a Monte Carlo process of generating multiple 
realizations of three-dimensional maps of K using sequential simulation methods (Deutsch and 
Journel, 1998). The maps are produced within the original domain, using the same mesh 
discretization, spatial statistics, and conditioning data as the original model, thus preserving the 
spatial distribution of the hydrogeologic categories and K at the same scale. The maps are then 
used to generate distributions of hydraulic head under hydraulic boundary conditions that 
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incorporate the uncertainty in head identified in the DDA for these boundaries. This represents 
an improvement over the 1999 model, which did not incorporate uncertainty in head boundary 
conditions. 

2.2.2. Model Domain 


The model is reduced in size from the original model because previous modeling 
demonstrated that areas east, west, and south of the Faultless land withdrawal area (UC-1) do not 
contribute to predictions of radionuclide migration. The new domain is 3.6 km long on each side 
(approximately twice the length of the land withdrawal area) and is centered over UC-1, and like 
the 1999 model, is aligned in the north-south direction (Figure 2.1). The domain covers the same 
1,350-m vertical section included in the 1999 model. Also consistent with the 1999 model, each 
cell in the uniform mesh is a cube having edge dimensions of 50 m. 
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Figure 2.1. Map view of model domain showing the finite difference mesh, the locations of the UC-1 
land withdrawal boundary and the related CNTA boreholes. 

2.2.3. Description of Hydrogeologic Heterogeneity 

Maps describing the three-dimensional spatial geometry of the hydrogeologic categories are 
produced using the sequential indicator simulation (SIS) algorithm (Alabert, 1987). The SIS 
simulation process, parameter values, and conditioning data are identical to those used in the 
1999 model, so the resulting indicator maps match the earlier results. The maps are generated to 
the full size of the 1999 model domain to incorporate all available indicator conditioning data 
(, i.e ., all the wells with lithologic and/or resistivity data) necessary to reproduce the large-scale 
hydrogeologic structure. The indicator maps are then trimmed to the smaller size of the new 
domain for use as conditioning for the K maps. 


4 































































































































Five hundred K maps are simulated for the new domain using the sequential Gaussian 
simulation (SGS) algorithm (Deutsch and Journel, 1998) in the same manner as was used in the 
1999 model, though the conditioning data are restricted to wells inside the new, smaller domain. 
Conditioning data that were utilized in the 1999 model, but lie outside the domain of the new 
model, do not affect the transport simulations because their distance exceeds the spatial 
correlation of K and they are not located in the direction of groundwater flow from the Faultless 
working point. All the parameter values are identical to those in the 1999 model with the 
exception of the vertical correlation lengths of K for the three hydrogeologic categories. In the 
1999 model, correlation lengths of 150 m were determined from spatial analysis of the CNTA 
packer tests, but values of 50 m were used in the base-case flow model to improve calibration. 
Reconfiguration of the head boundary conditions as uncertain and the utilization of a different 
calibration procedure in the current model (both described in subsequent sections of this report) 
allowed a return to the original values of 150 m. 

2.2.4. Modeling Codes 

The SISIM and SGSIM codes (Deutsch and Journel, 1998) are used for the sequential 
indicator and sequential Gaussian simulations, respectively. The source codes are unmodified, 
though the original random number generator has been replaced by a lagged Fibonnaci generator 
(LFG) (Knuth, 1981; Marsaglia, 1985). This generator is initialized for each Monte Carlo 
realization by setting the seed to the Monte Carlo realization number and is more effective than 
previously utilized methods for producing streams of independent random deviates. Other codes 
used in the model that require random number generators have also been adapted to use the LFG. 

The modular groundwater flow code MODFLOW-88 (McDonald and Harbough, 1988) is 
implemented to solve for the hydraulic head field using the finite difference approximation of the 
general flow equation 
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where K xx , K m and K zz , are values of K along the x, y, and z axes, respectively, (m/day), h is the 
potentiometric head (m), W is the volumetric flux per unit volume representing sources and/or 
sinks of water (m' 1 ), S s is the specific storage of the porous media (m 1 ), and t is time (day). 
Equation (2.1) is the most general fonn of the groundwater flow equation; in the current model, 
K is heterogeneous (as defined by the K maps), but is isotropic such that K = K xx = K yy = K zz . In 
addition, there are no sources or si nk s (/.<?., W = 0) and the flow problem is treated as steady state 
(i i.e ., dh/dt = 0). 

Three modules are used for the application of MODFLOW to Faultless. The Basic (BAS) 
module specifies the other modules that are to be used, the geometry of the model domain, the 
boundary conditions, and the time steps. The Block-Centered Flow (BCF) module handles the 
grid discretization, aquifer type, and hydraulic parameters required to solve the finite difference 
equations. The Preconditioned Conjugate-Gradient (PCG) module uses modified incomplete 
Cholesky preconditioning to efficiently solve the matrix of finite difference equations (Hill, 
1990). The only modifications made to the original code involved input and output: new routines 
were added to read the three-dimensional K maps and to save the maps of hydraulic head. 
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2.2.5. Hydraulic Boundary Conditions 

Other than the size of the domain, the groundwater flow model is configured in much the 
same way as the 1999 model. The northern, southern and bottom faces are specified head 
boundaries, and the east and west faces are no-flow boundaries (this configuration is based on 
the predominantly north-south flow patterns simulated in the 1999 model). Heads on the 
specified-head boundaries are treated as uncertain and are varied using a method adapted from 
the DDA modeling so that each realization is based on a unique, but equally likely, distribution 
of boundary heads. 

As was done for the 1999 model, boundary heads are estimated using head measurements 
made in the straddle-packed intervals of the CNTA exploratory wells UCe-20, UCe-18 and 
HTH-1. For the current model, however, the specified heads include uncertainty through 
incorporation of the variation observed in the multiple straddle-packed intervals. Just as in the 
1999 model, the process uses HTH-1 as the starting point, as it is located much closer to 
Faultless than the other wells. Estimates of head at the elevations of the top and base of the 
model at the location of HTH-1 were calculated in the 1999 model by vertical extrapolation of 
HTH-1 head values. The uncertainty in these values is quantified by the 95 percent confidence 
interval associated with the linear regression fit to the ten head measurements in HTH-1 (Pohll 
and Mihevc, 2000). Similarly, the horizontal gradients at the top and base of the model are 
determined using the uncertainty in the linear regression of head data from UCe-20, UCe-18 and 
HTH-1. For each Monte Carlo realization, heads at the top and bottom edge nodes of the 
northern and southern faces are extrapolated using these horizontal gradients and the heads 
estimated at the top and bottom at the location of HTH-1. The heads of the remaining nodes on 
the specified head boundaries are obtained by linear interpolation between the heads on the 
edges. 

Using this process, uncertainty in hydraulic head measurements translates into a large 
degree of variability in the hydraulic boundary conditions of the flow model. Consistent with the 
conceptual model, groundwater flow is directed downward at the northern end of the domain, 
upward at the southern end, and to the north along the base. 

2.2.6. Calibration of the Flow Model 


Calibration of the flow model is evaluated using the average of squared differences between 
the measured head h m and the simulated head h s at each of the ten straddle packer intervals in 
HTH-1. The root mean square error (RMSE) is calculated for each flow realization m using the 
expression 


RMSE = 


l/n 'Z( h m -K)- 


i =1 


0.5 


( 2 . 2 ) 


where n is the number of calibration targets. The RMSE ranges from 0.76 to 8.3 m, with a mean 
value of 1.7 m, for the full set of 500 Monte Carlo realizations. It should be noted that the 
calibration of the 1999 model was accomplished through adjustment of the vertical correlation 
length of K and the mean of K in the alluvium (Category 1) using measured hydraulic heads in 
the uppermost packer intervals of HTH-1, UCe-18 and UCe-20 as the calibration targets. No 
adjustments to the construction and parameterization are made in the current model with the 
exception of the value of the vertical correlation scale of K and the uncertainty incorporated in 
the head boundary conditions, as described previously. 
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Because the construction and parameterization of the current model are both essentially the 
same as the corresponding volume of the 1999 model, the range of simulated heads, flow 
directions and fluxes are similar. Comparison of the mean head profiles simulated by both 
models at the location of HTH-1 shows that the new model provides an improved match to the 
observed head distribution near UC-1 (Figure 2.2). Calibration of the 1999 model was designed 
to match heads over a larger region and as a result, simulated heads at HTFI-1 were not matched 
to the degree obtained in the current model. In addition, the closer proximity of the specified 
head boundaries in the current model results in increased control over heads simulated in the 
model interior. However, incorporation of uncertainty in the head boundary conditions results in 
a greater range of simulation results in the current model as compared to the 1999 model, as 
indicated by the wider range of the 95 percent confidence interval. 
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Figure 2.2. Profiles of hydraulic heads simulated by the 1999 model and the new model for well HTH-1, 
as compared to the measured values. 
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In a traditional stochastic numerical flow and transport model using Monte Carlo 
techniques, each of the realizations of flow receives equal weight in the transport analysis. 
However, it is clear from the range of simulated results that some of the realizations fit the field 
data better than others. In an effort to honor site-specific field information throughout the 
modeling process, those realizations that are in good agreement with the field data are given a 
greater relative weight in the transport modeling than those that are in poor agreement. 

The weighting procedure utilized here is the generalized likelihood uncertainty estimator 
(GLUE) (Beven and Binley, 1992) that extends Monte Carlo random sampling to incorporate the 
goodness-of-fit of each realization. The goodness-of-fit is quantified by the likelihood measure 

I(f|©) = |X(9 2 k (2.3) 

where 


a = (h —h ) 

V s m /1 


(2.4) 


and L(Y | 0) is the likelihood of the vector of outputs, Y given© is the vector of random inputs, 


h s is the simulated head at the point i, h m is the observed head at that point, and A is a likelihood 

shape factor. The choice of N is subjective though its value defines its relative function. As N 
approaches zero, the likelihood approaches unity and each simulation receives equal weight, as 
in the traditional Monte Carlo analysis. As N approaches infinity, the simulations with the lowest 
RMSE receive essentially all of the weight, which is analogous to an inverse solution. In this 
study, the value of N is assumed to be unity, which is a value typically used for this type of 
analysis (Beven and Binley, 1992; Freer et al., 1996; Morse et a/., in review). 


Each of the 500 flow realizations are weighted based on an application of the Bayes 
equation of the form: 


P(&\Y) = 


L(Y | 0)P(0) 
C 


(2.5) 


where P(&) is the prior probability of the input parameters produced by the Monte Carlo 
simulation, L(Y | 0) is the likelihood measure from Equation (2.3), C is a normalization constant, 

and P(Q \ Y ) is the posterior density. The posterior density is the probability of the input 
parameters occurring after taking into account the likelihood measure and is used to calculate the 
contaminant boundary as described in Section 4.0. Each of the 500 flow realizations receives a 
weight based on Equation (2.5), which is simply a normalization of the likelihood such that the 
sum of the weights for all realizations is unity. In this way, each transport realization is 
appropriately weighted according to the goodness-of-fit of its hydraulic head distribution to the 
available field measurements. 

2.3 Results of the Flow Model 

Table 2.1 presents summary statistics regarding simulated values of indicator category, K, 
heads at the four specified head edges, values of RMSE for simulated heads, and errors in the 
mass balance of the flow model. 
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Table 2.1. Uncertainty in selected parameters simulated by the flow model. 


Parameter 

Mean 

Standard Deviation 

PDF of Category 

Category 1 

0.23 

1.5xl0" 3 

Category 2 

0.48 

6.4xl0" 2 

Category 3 

0.29 

6.4xl0" 2 

Hydraulic Conductivity, log 10 [ K (m/d)] 

Category 1 

-0.72 

0.03 

Category 2 

-2.6 

0.006 

Category 3 

-0.46 

0.008 

Hydraulic Head on Boundaries, h s (m AMSL) 

North, Top 

1670.6 

5.9 

North, Base 

1652.4 

6.4 

South, Top 

1662.3 

2.3 

South, Base 

1668.6 

2.7 

RMSE (m) 

1.70 

0.74 

Mass Balance Error (In-Out, %) 

-0.001 

0.009 
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3.0 TRANSPORT MODEL 

3.1 Conceptualization of Radionuclide Transport 

Although no additional field work has been conducted at Faultless since the completion of 
the 1999 model, significant improvements have been made to the radionuclide transport code for 
simulating matrix diffusion processes and refined values of retardation parameters have been 
obtained through continued laboratory experiments. Furthermore, following the lead of the 
Faultless DDA model, many of the transport parameters that were treated as deterministic in the 
1999 model are now treated as uncertain. These parameters include effective porosity for each of 
the three hydrogeologic categories, porosity of the welded tuff matrix, matrix diffusion 
parameters, and dissolution from the nuclear melt glass. Additionally, sorption and geochemical 
release parameters are treated as uncertain for transport cases that include sorbing radionuclides. 

With the exceptions noted above, the general conceptualization of radionuclide transport at 
Faultless simulated in the 1999 model is maintained. A brief overview of the principal 
components of the conceptual transport model is included here; the reader is referred to 
Pohhnann et al. (1999) for complete descriptions. Aspects of the model that differ from the 
earlier model and the methods developed to incorporate them are presented in the Section 3.2. 

The contaminants of interest are the radionuclides produced by the Faultless test and the 
daughter products created through radioactive decay. All of the initial radionuclide mass is 
assumed to be located within the cavity; no prompt injection into the surrounding formations is 
included. The mass data for the radionuclides produced by Faultless remain classified and cannot 
be presented in a public document. Transport calculations are performed here using a unit value 
as the initial mass. The transport results based on unit values of initial mass are then converted to 
true mass in a classified companion document. To make the general results of the 1999 transport 
model publicly available, Pohhnann et al. (1999) estimated the initial mass of H using the 
highest concentration measured in water samples from post-test well UC-1-P-2SR and a generic 
relationship related to depth of burial (Glasstone and Dolan, 1977). The initial mass for H 
estimated in this way and used for the construction of the boundary maps in the current model is 
2.6 x 10 18 pCi. 

Radionuclides are apportioned according to their volatility among surface deposits and 
volume deposits in nuclear melt glass. Nuclides within the glass are released to the groundwater 
according to their glass dissolution behavior and melt glass characteristics while those forming 
surface deposits are immediately available for transport by groundwater. However, it is assumed 
that no migration occurs until the cavity and chimney have refilled with groundwater and 
hydraulic equilibrium conditions have been reestablished, which is estimated to occur 30 years 
after the nuclear test. Molecular diffusion is neglected during this infill period, presuming that 
the flow of groundwater toward the hydraulic sink is much larger than the migration of 
radionuclides away from the sink driven by concentration gradients (Pohhnann et al., 1999, 
Appendix 3). 

The radionuclides produced by the nuclear test are grouped into six transport classes based 
on their ratio of hydraulic release to geochemical release (glass dissolution) and their retardation 
factor as they are transported by advection through the hydrogeologic environment. Diffusion of 
radionuclides from fractures into surrounding matrix blocks in hydrogeologic Category 3 is also 
considered. Parameters for simulating retardation and diffusion processes are estimated from 
laboratory experiments using cores from CNTA exploratory boreholes. Nuclides are grouped 
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according to assumed general sorptive behavior and assigned the same retardation values 
because comprehensive radionuclide-specific data are not available for the Faultless 
environment. 

3.2 Transport Modeling Methods 

The random walk particle tracking (RWPT) method is used to simulate the transport and 
fate of radionuclides within the random velocity fields generated by the flow model. The RWPT 
code used in the 1999 model has been substantially updated to improve the handling of the 
spatial variability of dispersion and porosity, to incorporate a new approach for simulating the 
matrix diffusion process, and to include a number of other enhancements. 

Traditional random walk methods ( e.g ., Kinzelbach, 1988; Tompson and Gelhar, 1990) 
usually rely on the assumption that medium properties such as porosity, 9, and dispersion 
coefficient, D, are sufficiently smooth in space. Discontinuities in effective subsurface transport 
properties that may arise in discrete velocity fields of numerical groundwater flow models 
violate this smoothness assumption (LaBolle et al, 2000). Therefore, when 9 or D is 
discontinuous, these standard methods fail (LaBolle et al., 1996) because the gradient tenns of D 
and/or 9 cannot be formally defined. LaBolle et al. (2000) developed generalized stochastic 
differential equations applicable to the case of discontinuous coefficients (e.g., dispersion 
coefficients) and developed a new random walk method that numerically integrates these 
equations. That method is applicable for cases of abrupt changes in transport parameters and 
velocity values. The new random-walk equations proposed by LaBolle et al. (2000) can be 
written as 

=X, +[K(X f ,0]At + [2D(F(X / + 3X),t)Atf 2 Z (3 1} 

where the displacement vector 8X is defined as 

8X = [2D(K (X t ,t))AtJ^ 2 ■ Z (3.2) 

The approach evaluates the advective component of particle movement using the velocity at 
the current particle position, (x t , y t , z t ), and at time t. The dispersive component is performed 
using dispersive coefficients evaluated at an intermediate location, (x, + 8x, y, + by, z, + bz), 
where the increments 8x, by, and bz represent dispersive steps from the current location, (x h y t , z t ) 
to the intennediate location (x, + 8x, y t + by, z t + bz). The details of the approach and how it is 
incorporated in the RWPT code are described in Pohlmann et al. (2002). 

The transport approach described thus far is appropriate for a porous medium; however, the 
conceptualization of the flow system at Faultless includes highly fractured zones (Category 3, 
which has correspondingly high flow velocities) adjacent to unfractured porous zones 
(Categories 1 and 2). A continuum approach is applied in the sense that effective fracture 
properties (high K and low porosity) are randomly assigned to the model cells rather than direct 
incorporation of discrete fractures. As a result, particles are tracked through space in the same 
manner as for a porous medium, but they experience high flow velocities when they pass through 
a fracture cell. 

Matrix diffusion is incorporated using the physically based method presented by Liu et al. 
(2000) that uses particle transfer probabilities to describe the diffusion of particles into the matrix 
and back into the fractures. This method assumes that each grid cell contains two overlapping 
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continua: one represents the fractures containing the rapid migration of groundwater and solutes, 
the second represents the surrounding matrix blocks that have a larger quantity of water (owing 
to the higher porosity of the matrix), but where advection is negligible. The expression for the 
forward transition probability is given by Liu et al. (2000) as 

2 D* fm A 

Pfm ~ V f (\-A)(B-b) At ( 3 - 3 ) 

where Pf m is the probability that a particle in a fracture will diffuse into the matrix, D* fm is the 

hydrodynamic dispersion coefficient between fractures and matrix, A is the fracture surface area, 
Vj i s the volume of water in the fractures, /. is a shape factor expressing how the matrix volume 
increases as one approaches the center of the matrix block ( e.g., X = 1/3 for parallel fractures), B 
is the half spacing between the fractures, b is the fracture half-aperture, and At is the time step. 
Particles in fractures diffuse into the matrix according to Pf m where they are subjected to local 
advection within the matrix (if any; none is assumed here for Faultless), and dispersion/diffusion. 

Particles can also diffuse back from the matrix to the fractures according to the reverse 
transfer probability, P m f, where they are then moved at the local groundwater velocity within the 
fracture. This approach provides a finite limit on the size of matrix blocks, and thus limited, 
rather than infinite, diffusion capacity of the blocks. However, as presented by Liu et al. (2000), 
the transfer probabilities account only for nonreactive solute transport. To describe matrix 
diffusion of sorbing radionuclides, we have modified the method to incorporate equilibrium 
adsorption characterized by a linear isothenn, which transfonns the reverse transfer probability 
according to the appropriate retardation factor. Pohlmann et al. (2002) describe the modified 
approach and present the expression for the reverse transfer probability, which conservatively 
ignores sorption in the fractures, as 


mf ~ VJ1-A)(B-b) At - 4 ) 

where P m f is the probability that a particle in the matrix will diffuse back into a fracture, and V m 
is the volume of water in the matrix block. The selection of the appropriate transfer probabilities 
for the Faultless case is presented below in Section 3.3.4. 

3.3 Transport Parameters 

3.3.1 Radionuclide Release 

The rock, fission products, and components of the nuclear device that are vaporized by the 
tremendous heat and pressure of the nuclear explosion quickly begin to condense and coalesce 
into a melt glass that contains a significant portion of the radioactivity produced by the nuclear 
reaction. The slow dissolution of this glass by groundwater releases radionuclides to the flow 
system where they can be transported away from the cavity and chimney. The release process is 
described by the geochemical release rate, which is a function of the linear rate constant for mass 
transfer between the solid and liquid phases, the specific surface area of the glass matrix, and the 
gram formula weight of the glass matrix (Pohlmann et al., 1999). Uncertainty in the release rate 
is developed from a uniform logio distribution of specific surface area, with endpoints of 10 and 
100 cnr/gm (Bourcier et al., 2000), as was done in the DDA. 
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3.3.2 Porosity 

The ranges in effective porosity uncertainty for the three hydrogeologic categories are the 
same as those used in the DDA. The lognonnal distributions of effective porosity for the 
alluvium and tuffaceous sediments categories are derived from measurements of total porosity of 
CNTA cores, a reasonable approximation considering that these sediments are primarily 
nonindurated and nonwelded and thus incorporate flow and transport only within the matrix. On 
the other hand, flow and transport in the welded tuffs and rhyolites is assumed to occur primarily 
in fractures and faults and thus measurements of fracture porosity are required. The uncertainty 
in fracture porosity of Category 3 was obtained from the literature because the CNTA database 
does not include these data. This distribution is uniform logio and ranges between 0.001 and 
0 . 01 . 

In addition to the fracture porosity, the matrix porosity of the welded tuff is required to 
calculate two other transport parameters: the retardation in the matrix and the transfer probability 
for diffusion from the matrix to the fractures. The range in porosity of the Category 3 matrix is 
treated as a lognormal distribution based on the measurements of total porosity reported in the 
CNTA database for cores of welded tuff. Table 3.1 lists the statistical parameters used to 
describe the four porosity distributions. 


Table 3.1. Parameters used to describe the uncertainty in porosity for the three hydrogeologic categories 
and the welded tuff matrix. 


Category 

Distribution Type 

Mean 

Standard Deviation 

Alluvium 6 e 

Lognormal 

2.92 

0.30 

Tuffaceous Sediments 9 e 

Lognormal 

2.78 

0.42 

Welded Tuff 6 e 

Uniform logio 

-3.0 (minimum) 

-2.0 (maximum) 

Welded Tuff 6 m 

Lognormal 

2.68 

0.39 


3.3.3 Retardation 


Retardation of radionuclides through the chemical sorption process is approximated using a 
linear isotherm that is parameterized by a retardation coefficient in the transport model. This is 
the same bulk sorption approach that was applied in the 1999 and DDA models. The retardation 
coefficients are detennined for representative reactive solutes from laboratory batch experiments 
conducted at DRI using rock samples from each of the three hydrogeologic categories. The 
expression for the retardation coefficient is 


R = 1 + 


P h K d 

9 


(3.5) 


where p h is the bulk density of the sample, K d is the distribution coefficient determined as the 
ratio of the mass of sorbate sorbed per mass of sorbent to the aqueous concentration of sorbate, 
and 9 is the porosity. The values of bulk density are 2.04, 2.01, and 2.20 gm/cm for the 
alluvium, tuffaceous sediments, and welded tuff, respectively (Pohhnann et al., 1999). 

The K d values are taken as the slope of the linear regression of equilibrium concentration 

versus sorbed concentration for each solute at a pH of 8. These field-scale values oiK d are then 

adjusted to the scale of the model by considering the specific surface area of the adsorbent used 
in the experiments to that of the bulk aquifer materials as collected from the cores (as described 
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in Pohlmann et al., 1999). Uncertainty in K d is quantified through standard linear regression error 
analysis, which leads to a t-distribution for . The other random term in the retardation 

expression, porosity, is also varied through the uncertainty range appropriate for each 
hydrogeologic category to produce the final distributions of the retardation coefficient. Table 3.2 
lists the parameters used to define the retardation uncertainty in the alluvium, tuffaceous 
sediments and welded tuff matrix for the six classes of radionuclide transport. 


Table 3.2. Parameters used to describe the uncertainty in retardation for reactive solutes in the three 
hydrogeologic categories. The uncertainty ranges for effective porosity are listed in Table 3.3. 


Category 

Distribution Type 

Mean 

Standard Deviation 

Transport Class 3 

Alluvium K d 

h 

1.32xl0" 2 

1.30xl0" 3 

Tuffaceous Sediments K d 

ti 

1.16xl0" 3 

1.14xl0" 4 

Welded Tuff Matrix K d 

ti 

6.1 lxlO" 4 

2.22x1 O' 5 

Transport Class 4 

Alluvium K d 

ti 

4.87xl0" 5 

1.89xl0" 5 

Tuffaceous Sediments K d 

h 

1.98xl0" 5 

1.89xl0" 5 

Welded Tuff Matrix K d 

ti 

2.47xl0" 5 

1.89xl0" 5 

Transport Class 5 

Alluvium K d 

h 

3.27xl0" 4 

6.16xl0" 5 

Tuffaceous Sediments K d 

ti 

7.10xl0" 4 

6.49xl0" 5 

Welded Tuff Matrix K d 

ti 

3.13xl0" 4 

5.33xl0" 5 

Transport Class 6 

Alluvium K d 

ti 

6.20xl0" 3 

7.74xl0" 4 

Tuffaceous Sediments K d 

h 

3.03xl0" 3 

1.98xl0" 4 

Welded Tuff Matrix K d 

ti 

3.06xl0" 3 

1.51xl0" 4 


3.3.4 Matrix Diffusion 


As described in Section 3.2.1, particle transfer probabilities account for diffusion between 
the fractures and the rock matrix in the RWPT code. Unfortunately, the parameters required to 
compute the transfer probabilities expressed in Equations (3.3) and (3.4) are very difficult to 
estimate, and no field tracer test data are available at CNTA for calibration. In the absence of this 
information, we have estimated the transfer probabilities from the matrix diffusion parameter k 
using an empirical relationship between Pf m and k developed through numerical modeling studies 
by Hassan and Mohammed (2003). According to their work, the two parameters can be related 
through the expression 

P fm = 0.26/cAt (3.6) 


1 /9 . 

where k (days' ) is the matrix diffusion parameter and At is the time step used in the transport 
model. The matrix diffusion parameter is defined (Cvetkovic and Dagan, 1994) as 


k = 6, 


iPJ. 


(3.7) 


where 6 m is the matrix porosity, D m is the effective diffusion coefficient in the rock matrix 

(mVday), R m is the dimensionless retardation coefficient in the rock matrix, and b is the effective 
half aperture of the fractures (m). 


15 



Uncertainty in k was incorporated in the DDA through Monte Carlo sampling of Equation 
(3.7) with all of the included parameters assumed to be random. However, when the DDA values 
of k are substituted in Equation (3.6), the time steps must be on the order of 10' days or less to 
obtain values of Pf m that do not exceed the maximum allowable transfer probability of one. 
Considering that the total time for which the contaminant boundary maps are calculated is 1,000 
years, time steps of this length would result in RWPT simulations having prohibitively long 
computational times. 

An alternate and in fact more conservative approach used here is to generate the distribution 
for k using the mean value from the 1999 model and assign the appropriate range of uncertainty 
from the DDA model, followed by the generation of the random distribution in the forward 
transfer probability using Equation (3.6). Using field and laboratory data from Faultless where 
available and augmented by information from the literature, Pohhnann et al. (1999) estimated the 
deterministic value of k to be 0.03 for the nonsorbing radionuclide H. This value is used as the 
mean to generate a logio distribution of k, maintaining the relative variation about the mean 
defined by the standard deviation of the DDA k distribution by maintaining the coefficient of 
variation, cv, between distributions. The coefficient of variation is defined ascv = s/x , where s is 
the standard deviation and x is the mean of the distribution. The target values of the mean and 
standard deviation for the arithmetic k distribution are 0.03 and 0.04, respectively, and the 
corresponding mean and standard deviation of the logio distribution are adjusted to obtain these 
target values. 

The reverse transfer probabilities are obtained through combination of Equations (3.3) and 
(3.4) as 


p r. K, Mi 

U V, e,v, (3-8) 


where 0 m is the matrix porosity of the welded tuff, 6 f is the fracture porosity, and V t is the total 
volume of the matrix-fracture system. Equation (3.8) is solved for P /m to obtain 


0 f 

P = P — 

mf 1 fin q 


(3.9) 


The parameters describing the uncertainty in the matrix diffusion parameters are listed in Table 
3.3. 


Table 3.3. Uncertainty ranges for the transfer probabilities used to simulate matrix diffusion in the 
welded tuff hydrogeologic category. 


Parameter 

Distribution Type 

Mean 

Standard Deviation 

Forward Probability, P fm 

Lognormal 

0.23 

0.25 

Reverse Probability, P fm 

Lognormal 

0.009 

0.016 


Matrix diffusion is incorporated only in the fractured densely welded tuffs and rhyolites of 
Category 3; the other categories are considered unfractured, so matrix diffusion is not included. 
The spatial variability of matrix diffusion properties is not considered. When first released, the 
radionuclide source is assumed to exist only in the fractures. 
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3.4 Transport Model Configuration 

The Faultless transport calculations employ the random walk particle tracking method and 
the three-dimensional Darcy flux fields using the same grid discretization and domain size as the 
groundwater flow model. The radionuclide source is assumed to be the entire Faultless cavity, 
which is simulated in the model as a cube having edge lengths of 200 m, and within which 
particles are uniformly distributed. The values of the transport parameters that are not treated as 
uncertain are listed in Table 3.4. The time-step length for each realization is calculated within the 
RWPT code using the values of porosity for the three categories associated with that realization. 
Time-step lengths are chosen so that the Courant numbers for the realization are less than one to 
ensure that particles are not transported a distance equal to the dimension of one grid cell (50 m) 
in a single time step. 


Table 3.4. Values of parameters that describe the configuration of the transport model and values of 
deterministic parameters that are assigned the same value for all six transport classes. Time 
step length is variable because it is determined for each realization from the maximum 
velocity. 


Parameter 

Value 

Location of Source, Nevada Central Coordinates 

Easting (m) 

191,675 

Northing (m) 

431,075 

Elevation (m AMSL) 

885 

Edge Length of Source (m) 

200 

Source Mass, M 0 

1.0 

Infill Time (years) 

30 

Total Simulation Time (years) 

1,000 

Time Step Length, At (days) 

2.5 to 1,526.8 

Longitudinal Dispersivity 

0.05 

Transverse Dispersivity 

0.005 

Molecular Diffusion 

0.0 

Number of Realizations 

500 


The radionuclides chosen for inclusion in the model are grouped into six solute classes 
based on their ratio of hydraulic release to geochemical release, geochemical release coefficient, 
and retardation factor (Table 3.5). The transport of the radionuclides in each solute class is 
simulated as a group. Radioactive decay was not incorporated in the transport model, but was 
applied in post-processing for each individual nuclide. Each transport class is simulated with a 
unit source mass and then individual radionuclide responses are calculated in a classified 
environment where the true source mass is included in the analysis. The decay process, including 
in-growth, is included at this point for individual radionuclides. 


Table 3.5. Values of parameters specific to individual transport classes. 


Parameter 

Class 1 

Class 2 

Class 3 

Class 4 

Class 5 

Class 6 

Release Ratio 

Hydraulic Release/ 
Geochemical Release 
Number of Particles 
Radionuclides 

1 . 0/0.0 
20,000 

3 h, 14 c, 

85 Kr, 85 Rb 

0.5/0.5 
20,000 
36 C1, 129 I 

0.6/0.4 
20,000 

90 Sr, 90 Y, 
90 Zr 

0.05/0.95 

40,000 

"Tc 

0.8/0.2 
20,000 

137 Cs 

0.05/0.95 

40,000 

151 Sm, 152 Eu, 154 Eu, 
234 U, 238 U, 237 Np, 
239 Pu, 240 Pu, 241 Am 
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4.0 BOUNDARY DELINEATION 
4.1 Review of Radiation Dose Units 


Ionizing radiation is photon (including, high-frequency UV; X and gamma rays, and cosmic 
radiation) and alpha- and beta-particle emissions of sufficiently high energy ( e.g ., > 34 eV) that 
they interact with the atoms composing matter, stripping orbital electrons, and creating ion pairs. 
When ionizing radiation penetrates cells, it deposits energy and directly ionizes biological 
materials in this manner, and it can also create toxic agents from the ionization of cellular water 
molecules, such as free radicals and hydrogen peroxide. Even at low doses, radiation is 
considered to pose some likelihood of lifetime excess cancer risk ( i.e ., a stochastic effect), which 
can be quantitatively equated to either morbidity (defined as fatal and nonfatal cancers) or 
mortality (defined as fatal cancers only). 

The activity of a radioactive isotope {i.e., radionuclide) is defined as the number of atomic 
transformations per second. One atom transformed per second is equal to 1 Becquerel (Bq, in 
International System of Units, SI), and 3.7 x 10 10 Bq equates to 1 Curie (Ci, the special unit of 
activity used before SI units were adopted). However, a Bq of activity does not necessarily 
equate to the number of radiations emitted by the radioactive isotope as it decays. For example, 
tritium ( H) releases 1 [S particle per atom transformed (i.e.. 1 radiation/s) but cobalt-60 ( Co) 
releases 1 [1 particle and 2 y rays per atom transfonned (i.e., 3 radiations/s), and potassium-42 
( K) releases 1 p particle and 0.2 y rays per atom transfonned {i.e., 1.2 radiations/s). 

Because activity measures the quantity of atomic transformations per second for a 
radionuclide, and atomic transfonnations release ionizing radiation, it is the activity and not the 
mass of a radioactive isotope that is relevant for assessing health risk. Yet, modeling 
groundwater flow and radionuclide transport involves accounting for the mass of a radionuclide 
per unit volume of water (concentration). Specifically, unclassified transport simulations are 
performed using a unit mass. In this way, one can perform the majority of the calculations in an 
unclassified environment and then transform the results based on the classified source term in a 
classified environment. The mass of a radionuclide of concern per unit volume of water can be 
converted to its activity per unit volume of water {i.e., activity concentration) using the specific 
activity (SA) for that radionuclide. 

The specific activity {SA) is the relationship between mass {e.g., mole) and activity (Bq) for 
a radionuclide, such that SA = A/N [Bq/mole], where the decay constant (A.) is expressed as 
1/seconds, and equals 0.693/fi/2, and t m is the half life of the radionuclide (and conversely 
tin = 0.693/A.), which is the time for one-half of the parent atoms that are present for a 
radionuclide to decay to daughter products (progeny), and N = atoms/mole. Thus, the activity 
concentration for an individual radionuclide can be detennined from its activity using the 
relationship shown in Equation (4.1): 


Bq (atoms transformed/s) = 


Mole 

0.693 6.023 x l()~^ atoms 

-x - 

t , (s) mole 

1/2 


(4.1) 


Accordingly, knowing the normalized mass per unit volume in water {i.e., from the 
unclassified transport simulations), the half life, the source mass, and the atomic weight for a 
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radionuclide, the unclassified simulation results can be converted to classified activity 
concentrations for that radionuclide. Then, that value can be used to evaluate health effects from 
exposure to the radiation emitted by the radionuclides of concern in the water. 

The absorbed dose from exposure to ionizing radiation is expressed in units of gray (Gy, in 
SI units), where 1 Gy is defined as the energy per unit mass imparted to matter equal to 1 J/kg 
[and it is also equal to 100 rad, where rad is the International Commission on Radiation Units 
and Measurements (ICRU) special dimension used before SI units were adopted]. The dose 
equivalent is the product of the absorbed dose and a dimensionless radiation quality or weighting 
factor (Q, representing an adjustment factor for distinguishing between the biological 
effectiveness of different types of radiation; for example, Q = 1 for beta, gamma, and X-ray 
emissions, and Q = 20 for alpha emissions) and expressed in units of sievert (Sv), where 1 Sv 
represents the energy per unit mass imparted to biological material equal to 1 J/kg (and it is also 
equal to 100 rem, where rem is the ICRU special dimension used before SI units were adopted). 
The effective dose equivalent (ede) is a concept for converting internal exposures to an 
equivalent whole-body dose (by summing the products of tissue-weighting factors and 
equivalent doses to an organ) so that doses from internal exposures can be added to those from 
external exposures to determine a total dose (see Moeller, 1992). The ede also is expressed in 
units of Sv, and the tissue weighting factor represents the fractional contribution that an 
individual body organ can make, as a consequence of internal exposure to radiation, to the total 
lifetime excess cancer risk were the entire body ( i.e ., all organs) to be irradiated uniformly (EPA, 
1988). 

To account for dose that is being delivered to tissue(s) as long as a radionuclide is present 
inside the body and decaying, the committed effective dose equivalent (CEDE) was conceived. 
This derived value is also expressed in Sv, and is defined as the sum of all doses projected to be 
received in the future from the intake of radionuclides in the current year, and by convention, the 
period over which a dose is considered committed in a tissue is 50 y following intake, a value 
originally selected arbitrarily for application to occupational exposures (EPA, 1988). 

The CEDE agrees with modern metabolic and dosimetric understanding, and the CEDE per 
unit intake (Sv/Bq), or CEDE applicable “dose conversion factor,” for internal exposure by 
ingestion, has been calculated and published by EPA (1988) and can be used to establish intake 
and concentration guidance in a regulatory context. Such application is considered appropriate 
because for radionuclides with long effective half lives (i.e., long radioactive half lives in 
combination with a long residence time in the body, especially relative to a 50-y period of 
commitment or even a 70-y lifetime), the CEDE may overestimate the total dose that actually 
would be expected to occur over the 50-y commitment period (Moeller, 1992), or by implication, 
even a 70-y lifetime. 

In contrast with using the CEDE approach for deriving intake and concentration limits for 
radionuclides, there is another outdated method that is not as valid, but has regulatory precedent. 
For example, the National Primary Drinking Water Regulations for radionuclides (EPA 2000a) 
identify 0.04 mSv/y (4 mrem/y) as the total annual dose equivalent to an organ or the whole body 
that cannot be exceeded from internal exposure to beta particle and photon radioactivity, and 
then this annual dose equivalent, along with a 2 L/d consumption rate, is used as a basis for 
setting annual average activity-concentration limits for drinking water intake for beta-particle- 
and photon-emitting radionuclides (see EPA, 1976, for procedure). This method is best described 
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as the “critical-organ dose-limit approach” because it simply involves using an acceptable dose 
limit for a critical organ as a fundamental parameter for setting a concentration limit. 

Today, a far superior and more defensible procedure exists as an alternative to using the 
“critical-organ dose-limit approach,” or even applying CEDE dose conversion factors, for 
developing concentration limits for radionuclides. This new procedure capitalizes on the 
development of radionuclide-specific lifetime radiogenic cancer risk coefficients [expressed as 
either a cancer mortality (fatal only) risk per unit activity or as a cancer morbidity (fatal and 
nonfatal combined) risk per unit activity {i.e., mortality risk/Bq or morbidity risk/Bq)] for the 
U.S. population published by EPA (1999), and derived using “ ... state-of-the-art methods and 
models that take into account age and gender dependence of intake, metabolism, dosimetry, 
radiogenic risk, and competing causes of death in estimating the risks to health from internal or 
external exposure to radionuclides.” Furthermore, according to EPA (1999), the purpose of these 
cancer risk coefficients is intended “ ... to support rulemaking ...” and “ ... encouraged to 
promote consistency in risk assessment ...” especially by federal agencies. 

In fact, the published cancer risk coefficients are suited for use in prospective assessments 
of potential cancer risk related to long-term exposure to radionuclides in environmental media. 
Specifically, they have been tabulated for over 800 radionuclides individually (not categorically, 
as is done for MCLs, which also only address dose and not risk) and include values applicable to 
low-acute doses or low-dose rates from internal exposure through various media, including 
drinking water. Technically, these radionuclide-specific cancer-risk coefficients can be applied to 
estimate the lifetime excess mortality or morbidity cancer risk due to chronic exposure over a 
lifetime to a constant environmental activity concentration by an average individual in a 
stationary population in the U.S. Consequently, summing the products of drinking water activity 
concentrations predicted for specific radionuclides (Bq/L), a corresponding radionuclide-specific 
cancer-risk coefficient ( e.g ., to be exhaustive, cancer morbidity risk/Bq), and a conservative 
estimation of a lifetime exposure to drinking water {e.g., 2 L/d x 365 d/y x 70 y/lifetime), will 
yield a total lifetime excess cancer risk (in this example, cancer morbidity risk) for exposure to 
all of the radionuclides considered. 

The computed total lifetime excess cancer risk can then be compared to a de facto, “de 
minimus ” level of risk (i.e., so low it can be considered negligible) that can be selected from the 
lower range of risk identified by EPA as reasonable for establishing regulatory standards for 
drinking water contaminants [i.e., an excess lifetime cancer risk that does not exceed 10' 4 
(1/10,000) and ideally is less than 10' 6 (1/1,000,000)]. For example, limiting lifetime excess 
cancer morbidity risk to 10" (1 or less would be a legitimate unifying consideration with regard to 
establishing a reasonably well-conceived health-protective contaminant boundary perimeter for a 
CAU. 

Finally, a systematic quantitative uncertainty analysis, using either the full complex 
computational risk-coefficient model, or even a simpler abstraction, is a difficult task due to the 
large number of parameters and radionuclides (EPA, 1999). For this reason, the risk coefficients 
are taken as a group to be centrally located within their uncertainty range, even though this may 
not always be true for some individual radionuclide-specific cancer-risk coefficients (EPA, 

1999). Application of the published tabulated individual radionuclide-specific cancer-risk 
coefficients can be expected to confonn to accepted regulatory practice and produce results that 
would be considered sufficient. Therefore, this document will focus on the use of the published 
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radionuclide-specific cancer-risk coefficients without addressing the quantification of uncertainty 
related to the risk factors. 


4.2 Risk-based Analysis 

One metric for determining the extent of contaminant boundaries is morbidity cancer risk. 
Morbidity cancer risk is the most exhaustive endpoint, as it addresses both fatal and nonfatal 
cancers. EPA (2000a) identified an excess lifetime cancer risk that does not exceed 10‘ 4 
(1/10,000) and ideally less than 10‘ 6 (1/1,000,000) as a reasonable basis for establishing 
regulatory standards, and thus this range can represent a “de minimus ” level of risk (i.e., so low 
it can be considered negligible). For purposes of this document, limiting lifetime excess cancer 
morbidity risk to 10' 6 will be used as the limit of lifetime excess cancer morbidity risk for 
establishing the perimeter of a contaminant boundary for a CAU. The following strategy explains 
the procedure for calculating such a boundary. 

The parameters that define radionuclide-specific risk calculations are: 

1. Radionuclide activity concentration [c/xJJ] 

2. Ingestion exposure rate [ E ] 

3. Average life duration [ D ] 

4. The applicable radionuclide-specific morbidity cancer-risk coefficient [/y] 

The radionuclide activity concentration is derived from the flow and transport simulations 
and it is assumed that the majority of the uncertainty in the risk calculation is derived from this 
parameter. The ingestion exposure rate (2 L/day) and average life duration (70 years) are 
expected values as detennined by EPA (1999). The morbidity cancer-risk coefficients are 
derived by EPA (1999) and are shown in Table 4.1. It is important to note that risk coefficients 
are not available for argon-39, krypton-85 and gadolinium-150 and as such, cancer risk cannot be 
calculated for these radionuclides. 


The risk (R(x,t)) is calculated as: 

R(x,t)= 1 c{x,t\E]p\A (4.2) 

i=l 1 J 

where Cj(x,t) is the activity concentration for a radionuclide at a given position x and time t 
(pCi/L), E is the ingestion exposure rate (L/D), D is the average life duration (days), and r, is the 
risk coefficient for radionuclide i. 

Equation (4.2) is computed with the assumption that the activity concentration of a 
radionuclide in the environmental medium is constant over the entire lifetime of the exposed 
population. This assumption can be addressed in one of two ways. In aquifers with relatively low 
velocities, the concentrations can be collected from the transport model at 70-year intervals. 
Aquifers that exhibit rapid velocities require time averaging as: 



j + 70 


f \ 

x,t . 

v J J 


(4.3) 


22 



Table 4.1. The morbidity cancer-risk coefficients as derived by EPA (1999) and half-lives for 52 
radionuclides. Note that risk coefficients are not available for argon-39, krypton-85, and 
gadolium-150. 


Radionuclide 

Isotope Symbol 

Half life (t 1/2 ; y) 

Risk/pCi 

Tritium 

H-3 

1.23E+01 

5.07E-14 

Carbon-14 

C-14 

5.73E+03 

1.55E-12 

Sodium-22 

Na-22 

2.61E+00 

9.62E-12 

Aluminum-26 

Al-26 

7.30E+05 

1.73E-11 

Chlorine-36 

Cl-36 

3.01E+05 

3.30E-12 

Argon-39 

Ar-39 

2.69E+02 

n/a 

Potassium-40 

K-40 

1.28E+09 

2.47E-11 

Calcium-41 

Ca-41 

1.03E+05 

3.53E-13 

Manganese-54 

Mn-54 

8.55E-01 

2.28E-12 

Cobalt-60 

Co-60 

5.27E+00 

1.57E-11 

Nickel-59 

Ni-59 

7.60E+04 

2.74E-13 

Nickel-63 

Ni-63 

1.00E+02 

6.70E-13 

Krypton-85 

Kr-85 

1.07E+01 

n/a 

Strontium-90 

Sr-90 

2.91E+01 

5.59E-11 

Zirconium-93 

Zr-93 

1.50E+06 

1.11 E-12 

Niobium-93m 

Nb-93m 

1.61 E+01 

8.03E-13 

Niobium-94 

Nb-94 

2.00E+04 

7.77E-12 

Technetium-99 

Tc-99 

2.13E+05 

2.75E-12 

Ruthenium-106 

Ru-106 

1.02E+00 

4.22E-11 

Palladium-107 

Pd-107 

6.50E+06 

2.50E-13 

Cadmium-113m 

Cd-113m 

1.41E+01 

2.87E-11 

Tin-121m 

Sn-121m 

5.50E+01 

2.34E-12 

Tin-126 

Sn-126 

1.00E+05 

2.56E-11 

Antimony-125 

Sb-125 

2.76E+00 

4.37E-12 

Iodine-129 

1-129 

1.57E+07 

1.48E-10 

Cesium-134 

Cs-134 

2.07E+00 

4.22E-11 

Cesium-135 

Cs-135 

2.30E+06 

4.74E-12 

Cesium-137 

Cs-137 

3.02E+01 

3.04E-11 

Cerium-144 

Ce-144 

7.80E-01 

3.52E-11 

Samarium-151 

Sm-151 

9.00E+01 

5.55E-13 

Europium-150 

Eu-150 

3.60E+01 

4.33E-12 

Europium-152 

Eu-152 

1.35E+01 

6.07E-12 

Europium-154 

Eu-154 

8.59E+00 

1.03E-11 

Europium-155 

Eu-155 

4.71E+00 

1.90E-12 

Gadolinium-150 

Gd-150 

1.80E+06 

n/a 

Holmium-166m 

Hm-166m 

1.20E+03 

8.03E-12 

Thorium-232 

Th-232 

1.40E+10 

1.0 IE-10 

Uranium-232 

U-232 

7.00E+01 

2.92E-10 

Uranium-233 

U-233 

1.59E+05 

7.18E-11 

Uranium-234 

U-234 

2.46E+05 

7.07E-11 

Uranium-235 

U-235 

7.04E+08 

6.96E-11 

Uranium-236 

U-236 

2.34E+07 

6.70E-11 

Uranium-238 

U-238 

4.47E+09 

6.40E-11 

Neptunium-237 

Np-237 

2.14E+06 

6.18E-11 

Plutonium-238 

Pu-238 

8.77E+01 

1.31 E-10 

Plutonium-239 

Pu-239 

2.41E+04 

1.35E-10 

Plutonium-240 

Pu-240 

6.56E+03 

1.35E-10 

Plutonium-241 

Pu-241 

1.44E+01 

1.76E-12 

Plutonium-242 

Pu-242 

3.75E+05 

1.28E-10 

Americium-241 

Am-241 

4.33E+02 

1.04E-10 

Americium-243 

Am-243 

7.37E+03 

1.03E-10 

Curium-244 

Cm-244 

1.81E+01 

8.36E-11 
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where c -\ x ’ t ) is a running average of the activity concentration over a 70-yr period. Note that it 

may be necessary to run the simulations to t = 1,070 years to calculate the running average. In 
either case, the risk is calculated at 70-yr intervals at all locations x = (x,y,z). Thus, based upon 
the model solution c/(x,t) ,R(x,tj) can be computed at all discrete locations and at all times of 

the simulation. The locus of discrete points (or grid blocks) where R(x,tj) ever takes on a value 
greater than 10' 6 would then define a zone of contamination above the health-protective standard 
with respect to cancer morbidity risk. Conversely, the locus of discrete points where R(x,tj ) is 
equal to or less than 10' 6 would define a zone of contamination that is below the health- 
protective standard with respect to cancer morbidity risk. 

In application, one must account for the uncertainties in the flow and transport model, but 
Equation (4.2) does not explicitly account for uncertainty. The flow and transport simulations are 
performed in a Monte Carlo environment, which accounts for the parametric uncertainty inherent 
in all models. Thus, for each location x and selected time t, there will be a distribution of 
Cj(x,t) and R(x,t j ) values that need to be evaluated. The probabilistic interpretation (i.e., 95 
percent confidence intervals) is calculated to determine the contaminant boundary as described in 
the uncertainty analysis section below. 

4.3 Regulatory Based Analysis 

Another metric for determining the extent of contaminant boundaries is based on the SDWA 
standards. This is the metric specified in the FFACO. The regulatory standard is based on three 
categories of radionuclides: 

1. Alpha (a) emitters, excluding uranium 

2. Beta (P) and photon emitters 

3. Uranium 

If one or more of the categorical MCFs is exceeded in the groundwater, then the fluid 
parcel is considered to exceed the regulatory MCE. 

238 241 

For the alpha-emitting radionuclides (i.e., Pu, and Am), the sum of the activity 
concentrations of alpha-emitting radionuclides is calculated and compared to the 15 pCi/F 
standard. If the fluid parcel exceeds the 15 pCi/F standard, then the MCE for this fluid parcel is 
assumed to be violated. 

3 90 129 137 

For the P- and photon-emitting radionuclides, such as H (tritium), Sr, “ I, and Cs, a 
“sum-of-the-fractions” procedure is involved (see EPA, 2000c). First, simulated P- and photon- 
emitting radionuclides are detennined. A conversion table is then created (by procedures 
described in EPA, 1976), or consulted (see Huber, 2001, who applied the methodology described 
in EPA, 1976) to detennine the activity concentration corresponding to an annual dose 
equivalent of 4 mrem/y for each of the P- and photon-emitting radionuclides for which an 
activity concentration is predicted. Next, each radionuclide-specific p value ( e.g ., pCi/F) is 
divided by the activity concentration (e.g., pCi/F) equivalent to the annual dose equivalent of 
4 mrem/y for that particular radionuclide. This fraction represents the contribution of this 
radionuclide to the maximum allowable 4 mrem/y limit for all P and photon emitters present. 
These fractions are then summed for all beta and photon emitters for each fluid parcel (i.e., finite 
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difference cell within the model domain). If the sum exceeds unity, then the location is assumed 
to exceed the MCL. 

For uranium isotopes, the sum of mass concentrations of uranium isotopes (i.e., pg/L) 
present at each location is calculated. If the sum of mass concentrations exceeds the MCL for 
uranium (30 pg/L), then the MCL is violated. 

The categorical MCL evaluations are detennined at every location in the model domain at 
70-yr intervals, which represents the mean lifetime expectancy. The probabilistic interpretation 
(i.e., 95 percent confidence intervals) is calculated to detennine the contaminant boundary as 
described in the uncertainty analysis section below. 

4.4 Uncertainty Analysis 

If there were no uncertainty in the solute transport model, the risk-based or regulatory based 
analysis would provide explicit definition of the boundary maps. In reality though, the 
groundwater models contain a certain degree of uncertainty, which needs to be represented in an 
analysis of the boundary delineation. This was recognized in the FFACO and expressed as a 
requirement that the contaminant boundary be predicted at a 95 percent level of confidence. This 
requirement can actually be interpreted in two different, but complementary, ways to map the 
three-dimensional extent where groundwater has been contaminated and the associated 
uncertainty. One approach is to identify the region in which there is 95 percent certainty that 
contaminants exist and exceed the health risk or regulatory threshold. The boundary will then 
encircle the region that meets this criterion. The alternative approach is to identify the region 
where there is 95 percent certainty that contaminants do not pose a health risk and therefore 
encompasses the region that fails to meet the boundary criterion. In other words, the water that 
does not pose a health risk is external to the region enclosed by the boundary. 

Although the differences in calculation between these two approaches to computing a 
contaminant boundary may appear subtle, they are dramatically different in terms of their 
interpretation. The first approach provides an area in which there is a high degree of certainty 
that risk values are not lower than the calculated value, while the second approach provides an 
area in which there is a low probability that one will exceed a specified risk. 

These two types of boundaries are a direct result of the uncertainty analysis, but with 
different perspectives: contaminated zone versus uncontaminated zone. The approaches are 
statistical complements of one another: a boundary predicted by the model to have a 95 percent 
level of confidence in encompassing the contaminant is also the boundary where there is a 5 
percent assurance that the water external to it is clean, and a boundary encompassing an area that 
has a 5 percent confidence level that it contains the contaminant, has a 95 percent certainty that 
the area external is clean. The two approaches yield the same result if a 50 percent significance 
level is chosen. The description of the boundary in Appendix VI of the FFACO (as presented in 
the Introduction) is ambiguous as to which of these two perspectives should be used for the 
contaminant boundary. The approach presented here follows the second alternative, identifying 
the region where there is a 95 percent certainty that contaminants do not pose a health risk. 

The uncertainty in the groundwater flow and transport model is handled directly within a 
Monte Carlo-type approach. Monte Carlo methods describe the uncertainty in model predictions 
based upon an analysis of the uncertainty in model input parameters. Bayesian methods are used 
in combination with Monte Carlo techniques to weight each Monte Carlo realization based on its 
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ability to reproduce observed system behavior. The Bayesian methods (see Section 2.2.6) allow 
one to “calibrate” the stochastic model within the framework of the uncertainty analysis. 

The general approach to quantify the model uncertainty is described by the following steps: 

1. Determine all model parameters that are considered uncertain. It is important to note that 
uncertainty is derived from parameter measurement errors, spatial heterogeneity, and 
errors in model conceptualization each of which can be included in the uncertainty 
analysis. 

2. Quantify the distributions of all uncertain parameters. This is typically done via a 
specification of a probability distribution function (pdf) for each parameter. 

3. Perform numerical simulation of groundwater flow and transport with randomly sampled 
parameter distributions. 

4. Calculate the likelihood function ( e.g., weights) for each realization. According to 
Equation (2.5), the likelihood function is a measure of the “fit” between simulated and 
observed system behavior. Typically, hydraulic heads are used as the basis for calculating 
the likelihood measure. 

5. Quantify the pdf for the simulated solute transport. In the context of calculating boundary 
maps, a pdf of a risk-metric is created for each location in the model domain at 70-yr 
intervals. The pdf takes the form of probability exceedence of the 10‘ 6 risk threshold for 
the risk-based approach. 

6. Determine the 95 percent confidence threshold for the contaminant boundary. 

Steps 1 through 4 are described in sections 2.0 and 3.0 and the details of the pdf analysis and 
associated confidence thresholds ( i.e., steps 5 and 6) are provided below. 

It is important to note that the groundwater flow and transport simulations provide 
predictions of radionuclide transport under ambient conditions. The boundary maps are 
determined based upon an assumption that a human may become exposed to contaminated 
groundwater by drinking the groundwater at a rate of 2 L/day. Although a rate of 2 L/day is 
small, access to groundwater would typically be from a standard domestic well which, according 
to Nevada Law, can pump up to 1,800 gallons/day (-2,500 m /day). In most situations, domestic 
well pumping will not have a large impact on the rate and/or direction in which contaminants 
move, but under certain situations, this may not be the case. The appendix presents an analysis of 
the impact of well pumping on solute transport behavior. The delineation of the contaminant 
boundary may be dependent on the pumping regime and may require further investigation if 
someone wants to utilize the groundwater resources adjacent to the nuclear test. 

The quantification of the pdf for the risk-based analysis utilizes the results of the solute 
transport model and the source mass for individual radionuclides to detennine if a fluid parcel 
exceeds a 10" 6 risk threshold with the inclusion of model uncertainty. The risk-based approach 
yields a continuous pdf for each spatial position at 70-yr intervals. Although the statistical 
analysis presented herein is for a risk-based boundary, the analysis is similar for a regulatory 
based boundary. 

A confidence interval defines a range of risk values for a given significance level (e.g., 95 
percent), which describes the magnitude of the uncertainty inherent in the numerical model. A 
typical confidence region as defined by statisticians is given by: 
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(4.4) 


b 

J /O) = Y 

a 

where f(x) is the pdf, y is the significance level and a and b are the lower and upper bounds of 
the confidence region ( i.e., confidence levels). Figure 4.1 is a graphical representation of this 
type of definition. The dashed region represents an area equal to y. This type of definition 
captures the central portion of the pdf, which is not appropriate for the definition of boundary 
maps, as the tails represent very high risk and very low risk values. 



Figure 4.1. Probability distribution function showing how typical confidence levels are determined. 


For the purposes of contaminant boundaries, the choice is made to focus on the region of 
“clean” water and as such, the lower end (i.e., uncontaminated region) of the pdf is used in the 
determination of the confidence threshold: 


aj 

jffx) 


O 


= y 


(4.5) 


where aj is the threshold risk value that captures an area of size y. Figure 4.2 shows a graphical 
representation of the proposed contaminant boundary confidence threshold. The contaminant 
definition captures all transport realizations with risk values greater than aj and determines the 
threshold risk value that captures our knowledge to a level equal to y. If for example a y was 
chosen at a 95 percent level, then the threshold value represents a risk value in which there is 95 
percent confidence that the risk value is less than or equal to aj. 
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Figure 4.2. Probability distribution function showing how the risk threshold is calculated for the 
contaminant boundary. 


The boundary maps are created by detennining the threshold risk value for individual 
locations at 70-yr intervals.. The weighted risk values (i.e., the probabilities incorporate the 
Bayesian likelihoods) are used to create a pdf of risk. The threshold risk value is determined as 
shown in Figure 4.3. If the risk threshold is less than or equal to 10' 6 , or the MCL of any 
category is not exceeded at the 95 percent level, then the location is considered to be external to 
the contaminant boundary. This process is repeated for all locations within the model domain. 

All cells that did not meet the above criterion define the contaminant boundary region. 

Therefore, the contaminant boundary is a three-dimensional surface that encloses all cells that 
have a risk threshold greater than 10" 6 (or cells that violate the MCL). The boundary maps will be 
presented in three two-dimensional sections representing the x-y, x-z, and y-z Cartesian planes. 
The boundary sections are two-dimensional projections of the three-dimensional boundary 
surface. For example, the aerial view (x-y plane) represents the maximum extent of the three- 
dimensional boundary or volume as projected to the surface as shown in Figure 4.4. A location is 
considered to be within the boundary if at any x-y location there is at least one vertical position 
that is found to be within the boundary. Similar projections are made for the x-z and y-z planes. 

There are two options that are available to present the temporal history of the contaminant 
boundary. One option is to produce maps that represent an instantaneous view of the boundary 
for a given 70-yr time interval and the other option is to present the cumulative boundary from 
time zero to a given time period. The latter option is being proposed for the presentation of all 
boundary maps, as it represents the maximum extent of the boundary at any given time. 
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Figure 4.3. An example of the probability distribution function and how the risk threshold is calculated 
for the contaminant boundary. 
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Figure 4.4. An example of a three-dimensional contaminant boundary and the associated mapping to an 
x-y cross section. 
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5.0 HYPOTHETICAL EXAMPLE 


The methodologies presented above are applied to a hypothetical groundwater flow system 
for evaluation. The use of a hypothetical example allows one to test the methodology on a 
groundwater system in which one can ensure that transport is occurring. 

The hypothetical flow system is constructed as shown in Figure 5.1. A Gaussian hydraulic 
conductivity field is used with the mean adjusted such that conservative solutes move on average 
two-thirds of the distance from the source to the downgradient end of the model domain. Six 
transport classes were simulated to represent the different combinations of hydraulic/puddle 
glass release and geochemical sorption. The transport parameters used are shown in Table 5.1 
and are analogous to values used for other modeling studies of radionuclide transport (Pohhnann 
et aL, 1999). An unclassified source tenn mass as derived from the average of Area 19 and 20 
tests located on the Nevada Test Site is used for this hypothetical case (Smith, 2001). Thirty-nine 
(39) radionuclides are simulated because source masses were not available for the entire set. Risk 
calculations were performed on all 39 radionuclides, but only 31 were included in the regulatory 
based analysis, as MCLs are not available for the entire set (Table 5.2). A total of 100 
realizations were simulated to assess the impact of uncertainty of hydraulic conductivity spatial 
variability on the determination of contaminant boundary maps. Although other section maps can 
be calculated, only aerial boundary maps are presented for the hypothetical simulations. 


Constant Head Face 



Figure 5.1. Three-dimensional depiction of the hypothetical groundwater flow and transport model. 
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Table 5.1. Ratio of hydraulic versus glass release, retardation coefficients, and an example isotope as 
used in the hypothetical solute transport model. 


Class 

1 

2 

3 

4 

5 

6 

Hydraulic/Glass Release 

1. 0/0.0 

0.5/0.5 

0.6/0.4 

0.05/0.95 

0.8/0.2 

0.05/0.95 

Retardation 

1 

1 

10 

1 

10 

100 

Example 

3 H 

129j 

90 Sr 

" g Tc 

137 Cs 

234 U 


Table 5.2. Radionuclides, transport class, half life, source activity, risk factor, emission type, and MCL 
as used in the hypothetical transport model. 

Number Isotope Transport Class Half Life 

(ti/ 2 ; y) 

Source Activity 
(Ci) 

Risk Factor 
(Risk/pCi) 

Emission 

Type 

MCL 

l 

Tritium 

1 

1.23E+01 

9.19E+05 

5.07E-14 

beta 

20,000 pCi/L 

2 

Carbon-14 

1 

5.73E+03 

7.30E+00 

1.55E-12 

beta 

2,000 pCi/L 

3 

Aluminum-26 

1 

7.30E+05 

1.18E-04 

1.73E-11 

beta 

not available 

4 

Chlorine-36 

2 

3.01E+05 

2.81E+00 

3.30E-12 

beta 

700 pCi/L 

5 

Calcium-41 

2 

1.03E+05 

2.16E+01 

3.53E-13 

beta 

not available 

6 

Cobalt-60 

2 

5.27+00 

5.24E-01 

1.57E-11 

beta 

100 pCi/L 

7 

Nickel-59 

2 

7.60E+04 

5.54E+01 

2.74E-13 

beta 

300 pCi/L 

8 

Strontium-90 

3 

2.91E+01 

1.57E+04 

5.59E-11 

beta 

8 pCi/L 

9 

Zirconium-93 

3 

1.50E+06 

5.49E-01 

1.11E-12 

beta 

2,000 pCi/L 

10 

Niobium-93m 

3 

1.61E+01 

1.00E+02 

8.03E-13 

beta 

1,000 pCi/L 

11 

Niobium-94 

3 

2.00E+04 

2.28E+00 

7.77E-12 

beta 

not available 

12 

Technetium-99 

4 

2.13E+05 

4.03E+00 

2.75E-12 

beta 

900 pCi/L 

13 

Palladium-107 

4 

6.50E+06 

2.07E-02 

2.50E-13 

beta 

not available 

14 

Cadmium-113m 

4 

1.41E+01 

1.53E+01 

2.87E-11 

beta 

not available 

15 

Tin-121m 

4 

5.50E+01 

5.68E+01 

2.34E-12 

beta 

not available 

16 

Tin-126 

4 

1.00E+05 

6.46E-01 

2.56E-11 

beta 

not available 

17 

Iodine-129 

2 

1.57E+07 

1.24E-02 

1.48E-10 

beta 

1 pCi/L 

18 

Cesium-135 

4 

2.30E+06 

4.16E-01 

4.74E-12 

beta 

900 pCi/L 

19 

Cesium-137 

5 

3.02E+01 

1.99E+04 

3.04E-11 

beta 

200 pCi/L 

20 

Samarium-151 

6 

9.00E+01 

7.51E+02 

5.55E-13 

beta 

1,000 pCi/L 

21 

Europium-150 

6 

3.60E+01 

1.46E+01 

4.33E-12 

beta 

not available 

22 

Europium-152 

6 

1.35E+01 

4.32E+02 

6.07E-12 

beta 

200 pCi/L 

23 

Europium-154 

6 

8.59E+00 

2.04E+02 

1.03E-11 

beta 

60 pCi/L 

24 

Holmium-166m 

6 

1.20E+03 

5.89E-01 

8.03E-12 

beta 

90 pCi/L 

25 

Uranium-232 

6 

7.00E+01 

3.35E+00 

2.92E-10 

alpha & uranium 

30 mg/L 

26 

Uranium-233 

6 

1.59E+05 

2.25E+00 

7.18E-11 

alpha & uranium 

30 mg/L 

27 

Uranium-234 

6 

2.46E+05 

1.62E+00 

7.07E-11 

alpha & uranium 

30 mg/L 

28 

Uranium-235 

6 

7.04E+08 

2.18E-02 

6.96E-11 

alpha & uranium 

30 mg/L 

29 

Uranium-236 

6 

2.34E+07 

6.22E-02 

6.70E-11 

alpha & uranium 

30 mg/L 

30 

Uranium-238 

6 

4.47E+09 

2.89E-02 

6.40E-11 

alpha & uranium 

30 mg/L 

31 

Neptunium-237 

6 

2.14E+06 

4.81E-01 

6.18E-11 

alpha 

15 pCi/L 

32 

Plutonium-238 

6 

8.77E+01 

9.43E+01 

1.31 E-10 

alpha 

15 pCi/L 

33 

Plutonium-239 

6 

2.41E+04 

2.54E+02 

1.35E-10 

alpha 

15 pCi/L 

34 

Plutonium-240 

6 

6.56E+03 

8.16E+01 

1.35E-10 

alpha 

15 pCi/L 

35 

Plutonium-241 

6 

1.44E+01 

1.18E+03 

1.76E-12 

beta 

300 pCi/L 

36 

Plutonium-242 

6 

3.75E+05 

4.43E-02 

1.28E-10 

alpha 

15 pCi/L 

37 

Americium-241 

6 

4.33E+02 

6.14E+01 

1.04E-10 

alpha 

15 pCi/L 

38 

Americium-243 

6 

7.37E+03 

2.36E-03 

1.03E-10 

alpha 

15 pCi/L 

39 

Curium-244 

6 

1.81E+01 

3.92E+01 

8.36E-11 

alpha 

15 pCi/L 
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5.1 Regulatory Boundary 

The hypothetical model was used to calculate the regulatory based (MCL) contaminant 
boundary through a 1000-yr period. In an effort to better understand the contribution of the three 
regulatory categories (alpha, beta, and uranium) to the composite boundary, the individual 
category maps were calculated along with the composite boundary in Figure 5.2. As shown in 
Figure 5.2, the alpha and uranium boundary is restricted to the source area, and the beta 
boundary encompasses a much larger area. The composite boundary is dominated solely by the 
beta-emitting radionuclides. Upon further examination of specific radionuclides, it is clear that 
carbon-14 and iodine-129 contributed to the large areal extent of the beta boundary. Carbon-14 
has a large half life (5.73 x 10 3 years), is not assumed to be entrained in the puddle glass, and is 
assumed to be a non-reacting isotope. Iodine-129 is also conservative, has a large half life (1.57 
x 10 7 years) and although only 50 percent of its mass is assumed to be trapped in the puddle 
glass, the MCL (1 pCi/L) is the smallest of all of the beta-emitting radionuclides. 


Uranium Concentration > 30 ug/L 



0 10 20 30 40 

I Index 



10 20 30 40 


Alpha Concentration > 15 pCi/L 


Sum of Beta Ratio > 1.0 



I Index 



I Index 


Composite (alpha, beta and uranium) 


Figure 5.2. Regulatory contaminant boundary at the 95 percent significance level and 1,000 years for 
uranium, alpha, beta, and all (composite) categories as calculated from the hypothetical 
model. The boundaries are plotted for an x-y section with i and j corresponding to the x and y 
directions, respectively. 
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5.2 Risk-based Boundary 

The risk-based boundary was calculated through a 1,000-year period and compared to the 
regulatory based boundary in Figure 5.3. The two boundaries are identical for the hypothetical 
model using the source term as derived from the average of area 19 and 20 nuclear tests. As 
noted above, combinations of radionuclides in groundwater might yield estimates of lifetime 
excess cancer that could exceed a risk-based standard of 10" 6 , even though a categorical MCL 
may not be exceeded. In the hypothetical example, the magnitude of the source tenn is relatively 
large, thereby yielding large dissolved concentrations. The solute transport model has a limited 
ability to resolve concentration values, which is dependent on the source term magnitude and the 
number of particles used to simulate radionuclide migration. Therefore, the hypothetical model 
was unable to resolve concentrations down to the levels that define the MCL (~ 1 pCi/L). It 


Regulatory Based Boundary 


I Index 

Risk-based Boundary 


CD 

TD 


I Index 

Figure 5.3. Regulatory (top) and risk-based (bottom) contaminant boundaries at the 95 percent 
significance level and 1,000 years as calculated from the hypothetical model. The boundaries 
are plotted for an x-y section with i and j corresponding to the x and y directions, 
respectively. 
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stands to reason that if additional particles were used in the transport calculations, then minor 
differences between the regulatory and risk-based boundaries could be observed. 

5.3 Comparison of Statistical Approaches 

Figure 5.4 shows a comparison of the two possible statistical interpretations of the 
contaminant boundary for three significance levels (5, 50 and 95 percent). Because of the 
complementary nature of the two boundary types, the boundary excluding the area of clean water 
(proposed here for the contaminant boundary calculation) at a 95 significance level is identical to 
the boundary encompassing the contaminated water at 5 percent. As a direct result, the boundary 


Confidence Level = 95% 


0 10 20 30 40 

I Index 



10 20 30 40 



10 20 30 40 

Boundary encompassing 
contaminated water 


Confidence Level = 50% 


Confidence Level = 5% 


4CH 



10- 

0-1 - 1 -,-,-,- 

0 10 20 30 40 

Boundary excluding clean water 
(proposed contamination boundary) 


Figure 5.4. A comparison of the two statistical approaches at 5, 50 and 95 percent confidence levels. The 
boundaries are plotted for an x-y section with i and j corresponding to the x and y directions, 
respectively. The focus of the boundaries on the left side is identifying where the 
contaminated water exists, while the focus on the boundaries on the right side is focusing on 
where the water is uncontaminated. 
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developed to exclude the area of clean water will always yield a larger impacted region if the 
same significance level is used to calculate both approaches. Likewise, the boundaries are 
identical at the 50 percent significance level. 

5.4 Temporal Boundary Maps 

The boundary maps are calculated as an accumulation of risk through the time period of 
interest. Figure 5.5 shows the contaminant boundary as calculated at 70, 280, 560, and 1,000 
years. One could calculate an instantaneous boundary representing the risk at a specific time, but 
at least in the case of the hypothetical example, both the instantaneous and cumulative 
approaches yield the same result. The boundaries are equivalent because of the large number of 
radionuclides used to calculate the boundary. For example, if only a single radionuclide was used 
(< e.g tritium), then the instantaneous boundary area would increase and then decrease with time 
due to the relatively short half life. In this case, the two instantaneous and cumulative boundaries 
would be different, especially at greater times. When all of the radionuclides are used in the 
boundary calculation, there are numerous intersecting plumes due to sorption and glass 
dissolution processes. The plume intersections produce a temporal filtering effect that leads to 
similar instantaneous and cumulative boundaries. 



I Index I Index 



0 10 20 30 40 0 10 20 30 40 

I Index I Index 

Figure 5.5. A time sequence of the proposed contaminant boundary at the 95 percent significance level as 
calculated from the hypothetical groundwater flow and transport model. The boundaries are 
plotted for an x-y section with i and j corresponding to the x and y directions, respectively. 
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6.0 CONTAMINANT BOUNDARY FOR THE FAULTLESS UNDERGROUND 
NUCLEAR TEST, CAU 443 

The contaminant boundary for the Faultless underground nuclear test is determined using 
the revised groundwater flow and transport model presented in section 2. Both a risk-based and 
regulatory based calculation is presented. The methodology used for the calculation focuses on 
identifying the region where there is a 95 percent certainty that water is “clean” based on either 
the risk or MCL standard. This region is external (outside) to the boundary. The region internal 
to this boundary is 5 percent certain to contain water exceeding the risk or MCL standard. 

Contaminant boundaries using the risk-based calculation at a significance level of 95 percent are 
presented at 100 and 1,000 years after detonation in Figures 6.1 and 6.2, respectively. As 
described in Pohlmann et al. (1999), the calculated groundwater infill time is 30 years. Although 
other time periods could be presented, the contaminant boundary region does not change 
significantly during the 1,000-year simulation period. The contaminant boundary identifies the 
region where there is 95 percent certainty that contaminants do not pose a health risk external to 
the boundary delineation. 

Contaminant boundaries using the regulatory based calculation at a significance level of 95 
percent are presented at 100 and 1,000 years after detonation in Figures 6.3 and 6.4, respectively. 
Although the regulatory and risk-based calculations are based on slightly different risk 
thresholds, the boundary maps are identical. The similarities are likely due to the slow migration 
of radionuclides, which causes sharp contaminant fronts with large activity concentrations. In 
this case, the activity concentrations within the contaminant plume are well above the MCL 
standard and would therefore cause large risk values. 
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Figure 6.1. Contaminant boundary using the risk-based approach at a 95 percent significance level at 100 
years after detonation. 
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Figure 6.2. Contaminant boundary using the risk-based approach at a 95 percent significance level at 
1,000 years after detonation. 
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Figure 6.3. Contaminant boundary using the regulatory based approach at a 95 percent significance level 
at 100 years after detonation. 
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Figure 6.4. Contaminant boundary using the regulatory based approach at a 95 percent significance level 
at 1,000 years after detonation. 
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7.0 CONCLUSIONS 

This report documents the revisions to the CNTA groundwater flow and transport model, 
describes the methodology to calculate the contaminant boundary, and presents the boundary for 
the Faultless underground nuclear test. 

The contaminant boundary for the Faultless underground nuclear test at the Central Nevada 
Test Area is calculated using a newly developed groundwater flow and radionuclide transport 
model that incorporates aspects of both the original three-dimensional model (Pohlmann et al ., 
1999) and the two-dimensional model developed for the Faultless data decision analysis (DDA) 
(Pohll and Mihevc, 2000). The original three-dimensional model was approved by the State of 
Nevada conditional on additional uncertainty analysis, which was performed in the DDA. To be 
true to the state's intent for a more rigorous incorporation of uncertainty, DOE/Nevada 
Operations Office decided to enhance the original three-dimensional model with the uncertainty 
features developed in the DDA before proceeding with determining contaminant boundaries for 
Faultless. This new model includes the uncertainty in the three-dimensional spatial distribution 
of lithology and hydraulic conductivity from the 1999 model as well as the uncertainty in the 
other flow and transport parameters from the 2000 DDA model. Additionally, the new model 
focuses on a much smaller region than was included in the earlier models; that is, the subsurface 
within the UC-1 land withdrawal area where the 1999 model predicted radionuclide transport 
will occur over the next 1,000 years. 

These boundaries define over a 1,000-yr time period, with uncertainty quantified, the three- 
dimensional extent of groundwater that would be considered safe for domestic and municipal use 
around the nuclear test under ambient (non-pumped) conditions. Two metrics for defining “safe” 
are considered: a 10' 6 risk threshold, and a regulatory based threshold. In application to Faultless, 
there is no difference between the boundary calculated using risk as compared to the SDWA 
regulations. 


43 




8.0 REFERENCES 


Alabert, F., 1987. Stochastic Imaging of Spatial Distributions Using Hard and Soft Information. 
Master's Thesis, Stanford University, Stanford, CA. 

American Society of Mechanical Engineers Center for Research and Technology Development 
and Institute for Regulatory Science (ASME/IRS), 2001. Technical Peer Review Report, 
Report of the Review Panel: Strategy for Remediation of Groundwater Contamination at the 
Nevada Test Site, The American Society of Mechanical Engineers, Technical Publishing 
Department, New York, NY, CRTD-Vol. 62. 

Andricevic, R., J.I. Daniels, and R.L. Jacobson, 1994. “Radionuclide Migration Using a Travel 
Time Transport Approach and Its Application in Risk Analysis,” J. Hydrol., 163, 125-145. 

Beven, K.J., & A.M. Binley, 1992. The future of distributed models: Model calibration and 
uncertainty prediction. Hydrological Processes, 6:279-298. 

Bourcier, W., S. Roberts, D.K. Smith, S. Hulsey, L. Newton, A. Sawvel, C. Bruton, C. Papelis, 
W. Urn, C. Russell and J. Chapman, 2000. Determination of Reactive Surface Area of Melt 
Glass. Lawrence Livennore National Laboratory report UCRL-ID-145181, 86 p. 

Cvetkovic, V. and G. Dagan, 1994. Transport of kinetically sorbing solutes by steady random 
velocity in heterogeneous porous formations. J. FluidMech., 265, 189-215. 

Daniels, J.I., K.T. Bogen, and L.C. Hall, 2000. Analysis of Uncertainty and Variability in 
Exposure to Characterize Risk: Case Study Involving Trichlorethylene Groundwater 
Contamination at Beale Air Force Base in California. Water, Air, Soil Pollut., 123:273-298. 

Daniels, J.I., R. Andricevic, L.R. Anspaugh, and R.L. Jacobson, 1993. “4. Risk-Based Screening 
Analysis of Ground Water Contaminated by Radionuclides Introduced at the Nevada Test 
Site,” in Pilot Study Risk Assessment for Selected Problems at the Nevada Test Site (NTS), 

J.I. Daniels, Ed., Lawrence Livennore National Laboratory, Livennore, CA, 

UCRL-LR-113891. 

Deutsch, C.V. and A.G. Journel, 1998. GSLIB Geostatistical Software Library and User’s Guide, 
1998. Oxford University Press, 2nd edition. 

Lederal Lacility Agreement and Consent Order for Nevada (FLACO), 2000. Appendix VI. 
Corrective Action Strategy Section 3. Underground Test Area, in Nevada Federal Facility 
Agreement (FFACO), Agreement Between the U.S. Department of Energy (National Nuclear 
Security Administration), the State of Nevada (Nevada Division of Environmental 
Protection), and the U.S. Department of Defense (Defense Threat Reduction Agency), 
available from U.S. Department of Energy, National Nuclear Security Administration, 
Nevada Operations Office, Las Vegas, NV (December 7, 2000; Revision: 1); reissued on 
compact disc January 2002, to include revisions to any part of the FLACO that occurred 
since the July 2001 version was released. 

Freer, J. and K. Beven, 1996. Bayesian estimation of uncertainty in runoff prediction and the 
value of data: An application of the GLUE approach. Water Resour. Res., 32(7):2161-2173. 

Glasstone, S., & P.J. Dolan, 1977. The Effects of Nuclear Weapons (3rd ed.). U.S. Department 
of Defense and U.S. Department of Energy. 


45 



Hassan, A.E. and M.M. Mohammed, 2003. On using the random walk particle tracking method 
to simulate transport in single-continua porous media: a comparative study. Journal of 
Hydrology. 

Hill, M.C., 1990. Preconditioned Conjugate-Gradient 2 (PCG2), A computer program for solving 
ground-water flow equations. U.S. Geological Survey, Denver, Colorado. 

Huber, D., 2001. Derived Concentrations (pCi/L) of Beta and Photon Emitters in Drinking 

Water, Yielding a Dose of 4 mrem/y to theTotal Body or to Any Critical Organ as Defined in 
NBS Handbook 69, Standards and Risk Management Division, Office of Ground Water and 
Drinking Water, Washington, DC (Personnel E-mail Correspondence; Voice Telephone: 
202-564-4878/202-260-9566). 

Kinzelbach, W., 1988. The random walk method in pollutant transport simulation, In 

Groundwater Flow and Quality Modeling, eds. E. Custodio, A. Gurgui, J.P. Lobo Ferreira, D. 
Reidel, Norwell, Mass, 227-246. 

Knuth, D.E., 1981. The Art of Computer Programming Vol. 2: Seminumerical Methods (second 
edition), Addison-Wesley, Reading, Mass. 

LaBolle, E., G. Fogg and A.F.B. Tompson, 1996. Random-walk simulation of solute transport in 
heterogeneous porous media: Local mass-conservation problem and implementation 
methods. Water Resour. Res., 32:583-593. 

LaBolle, E., J. Quastel, G. Fogg and J. Gravner, 2000. Diffusion processes in a composite porous 
media and their integration by random walks: Generalized stochastic differential equations 
with discontinuous coefficients. Water Resour. Res., 36:651-662. 

Liu, H.H., G.S. Bodvarsson and L. Pan, 2000. Determination of particle transfer in random walk 
particle methods for fractured porous media. Water Resour. Res., 36(3):707-713. 

Marsaglia, G.A., 1985. A current view of random number generators, in Computational Science 
and Statistics: The Interface, ed. L. Balliard, Elsevier, Amsterdam. 

Maxwell, R.M., S.D. Pehnulder, A.F.B. Tompson, and W.E. Kastenberg, 1998. On the 
Development of a New Methodology for Groundwater-Driven Health Risk Assessment. 
Water Resour. Res., 34: 833-847. 

McDonald, M.G. and A.W. Harbough, 1988. A Modular Three-Dimensional Finite-Difference 
Ground-Water Flow Model. Techniques of Water Resources Investigations, Book 6, Chap. 
Al, U.S. Geological Survey. 

Moeller, D.W., 1992. “Chapter 10. Standards,” in Environmental Health (Harvard University 
Press, Cambridge, MA), pp. 225-228. 

Morse, B.S., G. Pohll, J. Huntington, and R. R. Castillo (In Press). Stochastic capture zone 
analysis of an arsenic-contaminated well using the GLUE methodology. Water Resource 
Research. 

Pohll, G., and T. Mihevc, 2000. Data Decision Analysis: Central Nevada Test Area. Desert 
Research Institute, Division of Hydrologic Sciences Report 45179. 

Pohll, G., J. Tracy, and F. Forsgren, 1998. Data Decision Analysis: Project Shoal. Desert 
Research Institute, Division of Hydrologic Sciences Report 45166. 


46 



Pohlmann, K., J. Chapman, A. Hassan and C. Papelis, 1999. Evaluation of Groundwater Flow 
and Transport at the Faultless Underground Nuclear Test, Central Nevada Test Area. 
Publication No. 45165, Desert Research Institute, Division of Hydrologic Sciences. 

Smith, David K., 2001. Unclassified Radiologic Source Tenn for Nevada Test Site Areas 19 and 
20. Fawrence Fivermore National Faboratory report UCRF-ID-141706, 4p. 

Tompson, A.F.B. and F.W. Gelhar, 1990. Numerical simulation of solute transport in three- 
dimensional randomly heterogeneous porous media. Water Resour. Res., 26(10):2541 -2562. 

United States Department of Energy, Nevada Operations Office, 2000. United States Nuclear 
Tests, July 1945 through September 1992, Fas Vegas, NV, DOE/NV-209-REV 15, 
December 2000, United States Department of Energy, Nevada Operations Office. 

United States Environmental Protection Agency (EPA), 1976. National Interim Primary 
Drinking Water Regulations, U.S. Environmental Protection Agency, Office of Water 
supply, Washington, DC, EPA-570-9-76-003; Reproduced by U.S. Department of Commerce 
National Technical Infonnation Service, Springfield, VA; NTIS Accession No. PB267630 
(1976); see Appendix B-Radionuclides (pp. 129-142) and Appendix IV (in Appendix B): 
Dosimetric Calculations for Man-Made Radioactivity, (pp. 152-159). 

United States Environmental Protection Agency (EPA), 1988. Federal Guidance Report No. 11: 
Limiting Values of Radionuclide Intake and Air Concentration and Dose Conversion Factors 
for Inhalation, Submersion, and Ingestion, Office of Radiation Programs, United States 
Environmental Protection Agency, Washington, DC, 

EPA-520/1 -88-020 (September 1988). 

United States Environmental Protection Agency (EPA), 1999. Federal Guidance Report No. 13: 
Cancer Risk Coefficients for Environmental Exposure to Radionuclides, Office of Radiation 
and Indoor Air, United States Environmental Protection Agency, Washington, DC, EPA 402- 
R-99-001. 

United States Environmental Protection Agency (EPA), 2000a. “Part II. Environmental 
Protection Agency, 40 CFR Parts 9, 141, and 142: National Primary Drinking Water 
Regulations; Radionuclides; Final Rule,” Fed. Regist. 65(236), 76708-76753 (December 7, 
2000 ). 

United States Environmental Protection Agency (EPA), 2000b. “Part IV. Environmental 
Protection Agency, 40 CFR Parts 9, 141, and 142: National Primary Drinking Water 
Regulations; Radionuclides; Notice of Data Availability; Proposed Rule,” Fed. Regist. 

65(78), 21576-21628 (April 21, 2000), see Table II-3, p. 21605-21614. 

United States Environmental Protection Agency (EPA), 2002. “Section II. SDWIS Reporting, 
Violation, Determination, and SNC Definitions,” in Implementation Guidance for 
Radionuclides, Washington, DC, EPA816-F-00-002. 


47 



APPENDIX: ASSESSMENT OF THE DEPENDENCE OF SOLUTE TRANSPORT 
MIGRATION ON DOMESTIC WELL PUMPING 


INTRODUCTION 

A hypothetical groundwater flow and transport model was constructed to illustrate the 
impact of domestic well pumping on solute transport migration. A steady-state, two-dimensional 
flow model was used with constant head boundary conditions on the left- and right-most model 
boundaries to induce a hydraulic head gradient from left to right (Figure A-l). The other two 
edges of the domain were set as no-flow boundaries. A Gaussian hydraulic conductivity field 
was constructed with the mean set such that conservative solutes would move across the domain 
in a 1,000-year period and have a log-variance of 1.0. Only one realization was perfonned, as the 
purpose was to test the impact of well pumping and not to address uncertainty issues. Four cases 
were simulated to represent conditions of large and small effective porosity and hydraulic 
conductivity and pumping versus no pumping conditions. The effective porosity and hydraulic 
conductivity was three orders of magnitude larger for simulations 1 and 2 as compared to 
simulations 3 and 4. Arbitrary length units were used for all simulations. The pumping rate for 
simulations 2 and 4 that included a well was set to 2 L/day, which would be much smaller than 
expected from a typical domestic well, but it illustrates the potential impacts on solute migration 
at the low end of the spectrum. 

The source was located approximately near the left edge of the model domain with the 
primary flow direction being from left to right (Figure A-l). The source was represented as an 
instantaneous release of a conservative solute with a unit concentration. Time histories of solute 
breakthrough to the well (either pumping or not) and the spatial distribution of solute 
concentration at various times within the 1,000-year simulation are presented to determine the 
importance of including local pumping within the flow simulation. 

RESULTS 

The spatial distribution of solute concentration is shown in Figure A2 at 250, 500, and 750 
years for simulations 1 and 2. As can be seen in Figure A2, the impact of the well is minimal 
when the effective porosity and hydraulic conductivity are large. Figure A3 shows the solute 
concentrations at 25, 90 and 260 years for simulations 3 and 4. Clearly, under conditions of small 
hydraulic conductivity and effective porosity, even small pumping rates have a large impact on 
solute transport migration. When the porosity is small, the area of influence from a pumping well 
becomes much larger. In simulation 4, the well is able to capture the entire solute plume. 

The impact of pumping on solute breakthrough is shown in Figure A4, which shows the 
breakthrough curves at the well location for all four simulations. In Figure A4(a), when the 
effective porosity is large, there are only minor differences in the breakthrough curve behavior. 
When the effective porosity is increased, the breakthrough is much more rapid and 
concentrations smaller when a pumping well is included (Figure A4(b)). The rapid breakthrough 
is attributed to increased hydraulic gradients induced by the well, which cause increased 
velocities. The smaller concentrations are attributed to the large capture zone, which causes 
dilution at the wellbore. 
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CONCLUSIONS 


The results of the hypothetical simulations show the potential for wells to cause significant 
changes to the solute transport migration behavior. The differences between including a well and 
not including one become more pronounced as the hydraulic conductivity and effective porosity 
decrease. 

The results suggest that contaminant boundary maps that are produced with the assumption 
that domestic wells have only a minor influence on radionuclide behavior may need to be 
modified under certain situations. 


No-Flow Boundary 



No-Flow Boundary 


Figure Al. Hypothetical groundwater flow and transport components for the two-dimensional flow 
model. 
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Figure A2. Spatial distribution of solute concentration for simulations 1 and 2 in the appendix. Images a, 
b, and c represent the concentrations for the no-well simulation at 250, 500, and 750 years, 
respectively. Images d, e, and f represent the concentrations for the simulation that included 
at pumping well at 250, 500, and 750 years, respectively. 
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Figure A3. Spatial distribution of solute concentration for simulations 3 and 4 in the appendix. Images a, 
b, and c represent the concentrations for the no-well simulation at 25, 90, and 260 years, 
respectively. Images d, e, and f represent the concentrations for the simulation that included 
at pumping well at 25, 90, and 260 years, respectively. 
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Figure A4. Solute breakthrough curves for simulations 1 and 2 (top) and 3 and 4 (bottom). 
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